function coalign_stack, images, xshift=xshift, yshift=yshift, xinits=xinits, $
	yinits=yinits, xrange=xrange, yrange=yrange, im_shifted=im_shifted, $
	magnification=magn, verbose=verbose, help=help
;+
; NAME:
;	COALIGN_STACK
; PURPOSE:
;       Find (and optionally apply) the optimal (x,y) offsets to maximize the
;       correlation of multiple images relative to the first image
; EXPLANATION:
;       Applies CORREL_IMAGES and CORRMAT_ANALYZE to a series of images,
;		each paired with the first supplied image.  Returns the an array of
;		the optimal (x,y) offsets for each member of the stack relative
;       to the first.
;
; CALLING SEQUENCE:
;       Result = COALIGN_STACK(images
;                [, XSHIFT=, YSHIFT=, XINITS=, YINITS=, XRANGE=, YRANGE=,
;                IM_SHIFTED=, MAGNIFICATION=, /VERBOSE, /HELP] )
;
; INPUTS:
;       images = a 3-D array, a stack of 2-D images with IMAGES[*,*,0] being
;                  the reference frame against which the others will be shifted.
;
; OPTIONAL INPUT KEYWORDS:
;       XSHIFT = the + & - shift, in un-magnified pixels, to be applied in
;                 X direction, default=7.
;                 Ignored if MAGNIFICATION is set > 1.
;       YSHIFT = the Y direction + & - shifting, default=7.
;
;       XINITS = a vector containing the initial X pixel offsets of each
;                  image in the stack, relative to the first image.  Should
;                  contain one element for each image except the first
;                  (reference) image.
;       YINITS = a vector containing the initial Y pixel offsets.
;
;       XRANGE = a two-element vector, [xmin, xmax], restricting the X range
;                  to consider when comparing images
;       YRANGE = [ymin, ymax], restricting the Y range to consider
;
;       MAGNIFICATION = option changes the resolution of images during
;                         computation.  Set to an integer > 1 to increase
;                         resolution and calculate alignment to fractional
;                         pixels; in this case, XSHIFT and YSHIFT are
;                         automatically set to MAGNIFICATION+2.
;                         If 0 < MAGNIFICATION < 1 then resolution
;                         is reduced, with a reduction factor of the closest
;                         integer to (1./MAGNIFICATION).
;
;       /VERBOSE causes the progress and results of analysis to be printed.
;
;       /HELP prints the calling sequence for this routine and immediately
;            returns.	
;
; OUTPUTS:
;       Result is a 2xN array, where N is the number of images stacked in
;         IMAGES.  RESULT[0,*] gives the optimal x offsets for each image,
;         and RESULT[1,*] gives optimal y offsets.  (Results[*,0] = [0,0])
;
; OPTIONAL OUTPUT KEYWORD:
;       IM_SHIFTED = the original image stack with each image shifted by
;                 the optimal (x,y) pixel shift to match the first image.
;
; PROCEDURE:
;       Loop over all images except the first.  For each image, call
;       CORREL_IMAGES to calculate the correlation matrix between current
;       image and first image.  Pass results to CORMAT_ANALYZE to get
;		optimal offsets.  If IM_SHIFTED keyword is set, apply shift to
;       appropriate image in the stack.
;       Return the optimal offset array.
;
; MODIFICATION HISTORY:
;       Written June 2005 by Nathaniel Cunningham, Univ. of Colorado
;       
;-

On_error,2

;                       Check parameters.

if keyword_set(help) then begin
	print,'USAGE:  coalign_stack (image array OR fits file list, '
	print,'         x_inital_offsets, y_initial_offsets '
	print,'         [,im_shifted = im_shifted] [,/verbose]'
	print,'         [,magnification = magnification] [,/help])'
	return,-1
endif

verbose=keyword_set(verbose)
if n_elements (magn) eq 0 then magn = 1

; check if images provided
if size(images,/TNAME) eq "UNDEFINED" then begin
; get a list of files to process
images = dialog_pickfile(/multi,filter='*.fits', $
	title='FITS images to align')
if images[0] eq '' then return,-1
endif

; check if inarr is a string array
if size(images,/TNAME) eq "STRING" then begin
	imlist=images
	nims = n_elements(images)
	; get the size of the array from the first image
	img = mrdfits(images[0],0,h1,/unsigned)
	imdims = size(img,/dimensions)
	ncol = imdims[0]
	nrow = imdims[1]
	tempimg = fltarr(ncol,nrow,nims,/NOZERO)
	tempimg[0,0,0] = temporary(img)
	for fn=1,nims-1 do begin
	   tempimg[*,*,fn] = mrdfits(images[fn],0,h1,/unsigned)
	endfor
	images = temporary(tempimg)
endif

;get size of input array
ndims=size(images,/n_dimensions)
if ndims ne 3 then message, 'IMAGES must be 3-D stack of 2-D images'
dims=size(images,/dimensions)
ncol=dims[0]
nrow=dims[1]
nims=dims[2]
if nims lt 2 then message, 'IMAGES must be a 3-D array of at least 2 images'

;make sure lengths of xinits, yinits matches (3rd dimension of IMAGES) - 1
if n_elements(xinits) eq 0 then xinits=fltarr(nims-1)
if n_elements(yinits) eq 0 then yinits=fltarr(nims-1)
inidims=[n_elements(xinits),n_elements(yinits)]
wh=where(inidims eq [nims-1,nims-1],inicount)
if inicount ne 2 then message, $
	'XINITS and YINITS must each have one less element'+$
	' than the number of images stacked in IMAGES'
	
if n_elements(xshift) eq 0 then xshift=7
if n_elements(yshift) eq 0 then yshift=7

if n_elements(xrange) eq 2 then begin
	wxr=where(xrange ge 0 and xrange lt ncol,wxcount)
	if wxcount eq 2 and xrange[0] lt xrange[1] then begin
		xmin=fix(xrange[0])
		xmax=fix(xrange[1])
	endif else begin
		xmin=0
		xmax=ncol-1
	endelse
endif else begin
	xmin=0
	xmax=ncol-1
endelse
if n_elements(yrange) eq 2 then begin
	wyr=where(yrange ge 0 and yrange lt nrow,wycount)
	if wycount eq 2 and yrange[0] lt yrange[1] then begin
		ymin=fix(yrange[0])
		ymax=fix(yrange[1])
	endif else begin
		ymin=0
		ymax=nrow-1
	endelse
endif else begin
	ymin=0
	ymax=nrow-1
endelse

magnify = fix(magn) > 1
reduce = round(1./magn) > 1
if magnify gt 1 and verbose eq 1 then $
	print, 'Magnification on; XSHIFT = YSHIFT = '+strtrim(magnify+2,2)

if magnify le 1 then xyoptim=intarr(2,nims) else xyoptim=fltarr(2,nims)
xyoptim[*,0]=[0.,0.]
if arg_present(im_shifted) then im_shifted=images

im0=images[xmin:xmax,ymin:ymax,0]
im0nan=where(finite(im0,/nan),i0count)
if i0count ne 0 then begin
	print,'WARNING -- NaN values encountered in image 0'
	print,'    setting these locations equal to 0.0'
	im0[im0nan]=0.
endif
for jim=1,nims-1 do begin
	im1=images[xmin:xmax,ymin:ymax,jim]
	im1nan=where(finite(im1,/nan),i1count)
	if i1count ne 0 then begin
		print,'WARNING -- NaN values encountered in image '+strtrim(jim,2)
		print,'    setting these locations equal to 0.0'
		im1[im1nan]=0.
	endif
	if magnify ge 2 then begin
		matrix=correl_images(im0,im1, $
			xshift=xshift, yshift=yshift, xoffset_b=xinits[jim-1], $
			yoffset_b=yinits[jim-1], magnification=magnify, monitor=verbose)
		corrmat_analyze, matrix, xbest, ybest, xoff_init=xinits[jim-1], $
			yoff_init=yinits[jim-1],magnification=magnify,print=verbose
	endif else if reduce ge 2 then begin
		matrix=correl_images(im0,im1, $
			xshift=xshift, yshift=yshift, xoffset_b=xinits[jim-1], $
			yoffset_b=yinits[jim-1], reduction=reduce, monitor=verbose)
		corrmat_analyze, matrix, xbest, ybest, xoff_init=xinits[jim-1], $
			yoff_init=yinits[jim-1],reduction=reduce, print=verbose
	endif else begin
		matrix=correl_images(im0,im1, $
			xshift=xshift, yshift=yshift, xoffset_b=xinits[jim-1], $
			yoffset_b=yinits[jim-1], monitor=verbose)
		corrmat_analyze, matrix, xbest, ybest, xoff_init=xinits[jim-1], $
			yoff_init=yinits[jim-1], print=verbose
	endelse
	xyoptim[*,jim]=[xbest,ybest]
	if n_elements(im_shifted) ne 0 then begin
		if magnify le 1 then begin
			im_shifted[0,0,jim]=shift(images[*,*,jim],xbest,ybest)
		endif else begin
			xlarge=magnify*ncol
			ylarge=magnify*nrow
			xoff=fix(xbest*magnify)
			yoff=fix(ybest*magnify)
			expanded=shift(rebin(images[*,*,jim],xlarge,ylarge,/sample),$
				xoff,yoff)
			im_shifted[0,0,jim]=rebin(expanded,ncol,nrow)
		endelse
	endif
endfor

if size(inlist,/TNAME) ne "UNDEFINED" then images=imlist
return,xyoptim
end